function vmd_make_ic(fname1d, fbase3d)
%
%   vmd_make_ic(fname1d, fbase3d)
%
%   Reads file with 1D-density distribution. Creates 3D IC.
%
%   Domain information and parameters are edited in the file.
%
%---------------------------------------------------------

%-- system size

%jx=8; jy=1; jz=32; jv=256; jdir=2; Lx=130; Ly=Lx;

%jx=8; jy=8; jz=32; jv=256; jdir=8; Lx=130; Ly=Lx;

jx=48; jy=48; jz=32; jv=256; jdir=8; Lx=900; Ly=Lx;

%-- physics parameters

u0=1.0; 
uAngle= 1/jdir; 
epsilon=1.e-4;


%---------------------------------------------------------
%-- reading 1D density distribution

  fid = fopen(fname1d, 'rb');
  dfu0 = fread(fid, jv*jz, 'double');
  fclose(fid);
  dfu0 = reshape(dfu0, jv,jz)/jdir;

%-- create and write 3D density --

  dfu = zeros(jv,jz,jy,jx);
  for ix=1:jx
     for iy=1:jy
         dfu(:,:,iy,ix) = dfu0;
     end
  end


  %for idir=1:jdir
  %fname = [fbase3d, '.dfu.0000.', num2str(idir-1,'%02d')];

  fname = [fbase3d, '.dfu.0000.00'];
  fid = fopen(fname, 'wb');
  fwrite(fid, dfu, 'double');
  fclose(fid);

  %end



%-- backround velocities --------------------------------

  theta = 2*pi*(0:jdir-1)/jdir + pi*uAngle;
  ux0 = u0 * cos(theta);
  uy0 = u0 * sin(theta);

  dx = Lx/jx;
  dy = Ly/jy;
  x  = (0:jx-1)*dx - Lx/2;
  y  = (0:jy-1)*dy - Ly/2;

  kx = 2*pi/Lx;
  ky = 2*pi/Ly;

  indx = round(jx/2-jx/6):round(jx/2+2+jx/6);
  indy = round(jy/2-jy/6):round(jy/2+2+jy/6);
  umod3d = zeros(jz,jy,jx);

  [xx,yy]= meshgrid(x, y);
  %noise = cos(xx*kx).*cos(yy*ky);

  %noise = cos(xx*kx);


  for idir=1:jdir


    noise = 2*rand(jy,jx)-1 + i*(2*rand(jy,jx)-1);
    noise(indy,:) = 0;  
    noise(:,indx) = 0;
    noise = real(fftn(noise))/sqrt(jx*jy)/2;


    umod2d = epsilon*noise + 1;
    for iz=1:jz
      umod3d(iz,:,:) = umod2d;
    end

    ux =  ux0(idir)*umod3d;

    fname = [fbase3d, '.ux.0000.', num2str(idir-1,'%02d')];
    fid = fopen(fname, 'wb');
    fwrite(fid, ux, 'double');
    fclose(fid);


    noise = 2*rand(jy,jx)-1 + i*(2*rand(jy,jx)-1);
    noise(indy,:) = 0;
    noise(:,indx) = 0;
    noise = real(fftn(noise))/sqrt(jx*jy)/2;

    umod2d = epsilon*noise + 1;
    for iz=1:jz
      umod3d(iz,:,:) = umod2d;
    end

    uy =  uy0(idir)*umod3d;
    fname = [fbase3d, '.uy.0000.', num2str(idir-1,'%02d')];
    fid = fopen(fname, 'wb');
    fwrite(fid, uy, 'double');
    fclose(fid);

  end


end






%---------------------------------------------------------
